Three practical abilities, in this order:
- Analytical benchmarks (queueing formulae) — the known answers against which you verify a simulator.
- Input generation — how to manufacture randomness of a required distribution from uniform numbers.
- Output interpretation — how much of the result to believe, and which inputs the result actually depends on.
Sections 1–3 are heavily numerical: expect derivations and small calculations in the examination, so practise the worked examples with a calculator.
1. Probability Refresher for Simulation
1.1 Random variables and their descriptors
A random variable X is described by its cumulative distribution function F(x ) = P(X ≤ x), and by its pmf p(x) (discrete) or pdf f(x) (continuous). The descriptors you will use constantly:
The coefficient of variation CV is the key number in queueing: it measures variability on a scale that is independent of units, and it is CV, not the mean, that determines how bad a queue becomes.
1.2 Distributions you must know by name and shape
| Distribution | Mean / variance | Typical simulation use |
|---|---|---|
| Uniform(a,b) | (a+b)/2 ; (b−a)2 /12 | Complete ignorance between bounds; the source of all other variates |
| Exponential(λ) | 1/λ ; 1/λ2 (CV = 1) | Inter-arrival times of a Poisson process; memoryless service; time to failure |
| Poisson(λt) | λt ; λt | Number of arrivals in an interval |
| Erlang(k,λ) / Gamma | k/λ ; k/λ2 (CV = 1/√k) | Service made of k sequential phases; less variable than exponential |
| Normal(μ,σ2) | μ ; σ2 | Measurement error, sums of many effects; never for a non-negative duration without truncation |
| Lognormal | eμ+σ2/2 | Task durations, file sizes, repair times — right-skewed and positive |
| Weibull(k,λ) | λΓ(1+1/k) | Reliability and lifetimes with increasing or decreasing hazard rate |
| Triangular(a,m,b) | (a+m+b)/3 | Expert opinion when only min, mode and max are known |
| Empirical | From the data itself | When you have plenty of measurements and no theoretical form |
Using Normal for a service time or an inter-arrival time. The Normal distribution extends to −∞, so it will eventually generate a negative duration, and your simulation will either crash or silently schedule an event in the past. Use lognormal, gamma or a truncated distribution instead.
1.3 Input modelling: fitting a distribution to data
- Plot the data — histogram, box plot, time series (to check for trends and non-stationarity).
- Choose candidate families from the shape and from physical reasoning.
- Estimate parameters (maximum likelihood, or method of moments).
- Assess the fit: Q–Q plot, P–P plot, and a goodness-of-fit test — chi-square, Kolmogorov–Smirnov, or Anderson–Darling (which weights the tails, usually the part that matters).
- If nothing fits, use the empirical distribution — but then you can never generate a value outside the observed range, which matters for rare events.
2. Introduction to Queues
A queueing system consists of an arrival process, one or more servers with a service-time distribution, a waiting room of finite or infinite capacity, and a queue discipline that decides which waiting customer is served next.
2.1 Kendall's notation
A queue is described as A/S/c/K/N/D:
- A — inter-arrival distribution: M (Markovian/exponential), D (deterministic), Ek (Erlang), G (general).
- S — service distribution, same codes.
- c — number of parallel servers.
- K — system capacity (default ∞), N — population size (default ∞), D — discipline (FCFS by default; also LCFS, SIRO, priority, processor sharing).
So M/M/1 is Poisson arrivals, exponential service, one server, infinite queue, FCFS; M/D/3/10 is Poisson arrivals, constant service, three servers and a waiting room of ten.
2.2 Performance measures and Little's law
Standard notation: λ = arrival rate, μ = service rate per server, ρ = λ/(c μ) = utilisation, L = mean number in system, Lq = mean number in queue, W = mean time in system, Wq = mean waiting time.
It holds for any stable queueing system in steady state — any arrival process, any service distribution, any number of servers, any discipline. That generality is what makes it so useful: it converts a measured count into a measured time, and it is the first sanity check to run on simulation output.
2.3 The M/M/1 queue
For ρ = λ/μ < 1 (the stability condition), the steady-state probability of n customers is Pn = (1−ρ)ρn, and:
Jobs arrive at λ = 12/hour; mean printing time is 4 minutes, so μ = 15/hour.
- ρ = 12/15 = 0.8 → the printer is busy 80% of the time.
- L = 0.8/0.2 = 4 jobs in the system on average.
- W = 1/(15−12) = 1/3 hour = 20 minutes total.
- Wq = 20 − 4 = 16 minutes waiting.
- Check with Little: L = λW = 12 × 1/3 = 4. ✓
Now raise the load to λ = 14: ρ = 0.933, L = 14, W = 60 minutes. A 17% increase in load tripled the queue. That non-linearity is the single most important practical lesson of queueing theory.
2.4 Other standard results
- M/M/c — multi-server; the waiting probability is given by the Erlang-C formula. Pooling servers is always better than separate queues: one queue with c servers beats c queues with one server each.
- M/M/c/c — no waiting room; blocking probability from the Erlang-B formula. The classical telephone-trunk model.
- M/G/1 — general service, given by the Pollaczek–Khinchine formula:
Read the second form carefully: with the same mean service time, halving the variability (CV from 1 to 0) halves the queue. Deterministic service (M/D/1) has exactly half the waiting of M/M/1. Variability, not just load, causes congestion — and this is precisely what a simulation captures and a spreadsheet of averages does not.
2.5 Queueing networks and why we simulate
Jackson networks (open, Poisson arrivals, exponential service, probabilistic routing) decompose into independent M/M/c nodes, which is elegant but demands assumptions that real systems violate: finite buffers with blocking, priorities, batch arrivals, correlated inter-arrivals, time-varying rates, balking and reneging, server breakdowns. Once any of these is present the closed form disappears and simulation is the only option. Use the analytical results as verification benchmarks: simulate M/M/1, confirm that the measured L and W match the formulae within the confidence interval, and only then trust the simulator on the harder model.
3. Random Noise
Noise is the random component of a signal or of a system's behaviour — the part not explained by the deterministic model. In simulation, noise is both something we inject (to model real disturbance) and something we suffer (as sampling error in the output).
3.1 Types of noise by spectrum
| Noise | Power spectral density | Where it appears |
|---|---|---|
| White | Flat, S(f) ∝ f0 | Thermal/Johnson noise, independent measurement error; samples uncorrelated |
| Pink (1/f) | S(f) ∝ 1/f | Flicker noise in electronics, network traffic, human-generated processes |
| Brown/red (random walk) | S(f) ∝ 1/f2 | Integrated white noise; drift in sensors, Brownian motion, asset prices |
| Shot | Poisson arrivals of discrete quanta | Photon or electron counting; packet arrivals |
| Impulse / burst | Rare, large deviations | Equipment faults, outliers, switching transients |
3.2 Modelling noise in a simulation
- Additive vs multiplicative: y = x + ε models a fixed-magnitude disturbance; y = x(1+ε) models a proportional one. Choosing wrongly distorts the tails badly.
- White Gaussian noise in a continuous model turns the ODE into the stochastic differential equation dX = a(X,t)dt + b(X,t)dW, integrated by Euler–Maruyama:
Note the √Δt: the noise term scales with the square root of the step, not with the step. Getting this wrong makes the noise magnitude depend on your step size — a classic bug.
- Coloured noise is generated by filtering white noise, most simply with an autoregressive process xt = φxt −1 + εt, which introduces the temporal correlation real disturbances have.
- Signal-to-noise ratio SNR = 10 log10(Psignal /Pnoise) dB is the standard way to specify how much noise to inject.
3.3 Pseudo-random number generation
A pseudo-random number generator (PRNG) is a deterministic algorithm xn+1 = g(xn) which, from a seed, produces a sequence un ∈ (0,1) that is statistically indistinguishable from independent uniform draws for the purposes of the simulation.
Required properties: long period (at least the square of the numbers you will consume), good uniformity and independence, high-dimensional equidistribution, speed, small state, reproducibility, and independent substreams for parallel runs (Unit III).
The classical linear congruential generator is xn+1 = (axn + c) mod m, with un = xn/m. It is easy to analyse and easy to get wrong: LCG points fall on a limited number of hyperplanes in d dimensions (Marsaglia's theorem), and the notorious RANDU generator produced only 15 planes in three dimensions, invalidating a generation of published results. Modern practice uses Mersenne Twister (period 219937−1), MRG32k3a (designed for streams), PCG or xoshiro. Testing is done with the TestU01, Dieharder or NIST batteries.
Never use a cryptographic-strength requirement or a language's default `rand()` for serious stochastic simulation without checking its period and quality; and never seed with the clock when you need reproducibility. Every reported result should be accompanied by the generator name, the seed and the substream assignment.
4. Random Variates Generation
Random variate generation is the process of transforming uniform (0,1) pseudo-random numbers into observations from a specified probability distribution.
4.1 Inverse transform method
If U ~ Uniform(0,1) and F is a cdf, then X = F−1(U) has cdf F.
Proof: P(X ≤ x) = P(F−1(U) ≤ x) = P(U ≤ F(x)) = F(x), because F is non-decreasing and U is uniform. ■
F(x) = 1 − e−λx. Set u = 1 − e−λx and solve:
(The second form is valid because 1−U is also Uniform(0,1); guard against U = 0.) With λ = 0.5 and u = 0.23, x = −2 ln(0.77) = 0.523.
Other inverse-transform cases you should be able to derive on demand:
- Uniform(a,b): X = a + (b−a)U.
- Weibull(k,λ): X = λ(−ln U)1/k.
- Triangular, Pareto, logistic — all have closed-form inverses.
- Discrete distributions: build the cumulative table and return the first xi with F(xi) ≥ u; use a binary search or the alias method for large supports.
- Empirical continuous: linear interpolation of the empirical cdf.
4.2 Acceptance–rejection method
Used when F−1 has no closed form. Choose a majorising function cg(x) ≥ f(x) where g is easy to sample:
- Generate Y ~ g.
- Generate U ~ Uniform(0,1).
- If U ≤ f(Y)/(cg(Y)) accept X = Y; otherwise reject and repeat.
The expected number of trials per accepted variate is exactly c, so the method is efficient only when the envelope hugs the density.
Take g = Uniform(0,1) so g(x) = 1, and c = 3 since max f = 3. Accept Y when U ≤ Y2. On average three pairs are drawn per accepted value. (Here the inverse transform X = U1/3 is obviously better — which is the point: use rejection only when inversion is unavailable.)
4.3 Special-purpose methods
- Convolution — a sum of variates: Erlang(k,λ) = −(1/λ)ln(U1…Uk); Binomial as a sum of Bernoullis.
- Composition — for a mixture distribution, first pick the component, then sample from it.
- Box–Muller for the Normal: Z1 = √(−2 ln U1) cos(2πU2 ), Z2 = √(−2 ln U1) sin(2πU2) — two independent standard normals per pair of uniforms; the polar (Marsaglia) variant avoids the trigonometric functions. Then X = μ + σZ.
- Poisson process — generate exponential gaps; for a non-homogeneous process with rate λ(t), use thinning: generate at the maximum rate λ* and accept each point with probability λ(t)/λ*.
- Correlated variates — use the Cholesky factor of the covariance matrix: X = μ + LZ, or copulas for non-Normal dependence.
4.4 Output analysis: making the numbers mean something
Generating the input is half of the statistics; interpreting the output is the other half.
- Terminating vs steady-state simulation. A terminating simulation has a natural end (a bank day, one mission); a steady-state simulation studies long-run behaviour and needs a warm-up (transient) period to be discarded — identified with Welch's graphical procedure.
- Independent replications. Run n replications with different substreams; the replication means are i.i.d., so the classical confidence interval applies:
- Batch means. For one long steady-state run, split the series into k large batches, treat the batch means as approximately independent, and apply the same formula. Batches must be long enough for the autocorrelation to decay — otherwise the interval is too narrow and you will be over-confident.
- Comparing two systems. Build a confidence interval on the difference. If it excludes zero, the difference is significant. Use common random numbers to induce positive correlation and shrink the interval.
- Choosing the run length. To reach a target half-width h* from a pilot run of n0 replications with half-width h0, use n ≈ n0(h0/h*)2 .
Computing a confidence interval from the individual observations within one run (e.g. the 10 000 customer waiting times). Those are strongly autocorrelated — a long wait is followed by long waits — so the usual formula badly underestimates the variance. Use replications or batch means. This error is common enough that examiners ask about it directly.
5. Sensitivity Analysis
Sensitivity analysis (SA) is the study of how the variation in a model's output can be apportioned to variation in its inputs — parameters, initial conditions, distributional assumptions and structural choices.
It is distinct from uncertainty analysis, which quantifies how much uncertainty the output has; SA asks where that uncertainty comes from.
5.1 Why it is done
- Factor prioritisation — which parameters deserve better measurement?
- Factor fixing — which parameters can be frozen at a nominal value, simplifying the model?
- Model verification — a parameter that should matter but does not (or vice versa) usually indicates a bug.
- Robustness of the recommendation — does the preferred policy remain the best across the plausible parameter range?
- Credibility — SA is part of validation and of any serious VV&A report (Units I and V).
5.2 Local (one-factor-at-a-time) sensitivity
Perturb one input around its nominal value and measure the response:
The relative (elasticity) form is what you should report, because it is dimensionless and comparable across parameters: “a 1% increase in service rate reduces mean wait by 4.6%”.
Limitations: local SA explores a single point of the input space, assumes linearity, and completely misses interactions between factors. For a stochastic simulation, Δ must also be large enough that the response rises above the noise — use common random numbers so that the difference is not swamped by sampling variation.
5.3 Global sensitivity analysis
Global methods vary all inputs simultaneously over their full ranges.
- Scatter plots and regression. Sample the inputs (Latin hypercube), run the model, and regress the output on standardised inputs. The standardised regression coefficients and the partial rank correlation coefficients (PRCC) rank the factors; valid when the model is monotonic and the regression R2 is high.
- Morris elementary effects (screening). A cheap one-step-at-a-time design over a grid, repeated from many random starting points. Reports μ* (mean absolute effect, = importance) and σ (spread, = non-linearity or interaction). Ideal as a first pass when the model has 20–100 inputs.
- Variance-based (Sobol) indices. The gold standard. Decompose the output variance into contributions from each factor and each interaction:
Si is the first-order index (the effect of xi alone); STi is the total index (its effect including all interactions). Interpretation rules worth memorising:
- ∑Si = 1 means the model is additive; a shortfall measures the strength of interactions.
- STi ≈ 0 means factor i can be safely fixed.
- STi » Si means factor i acts mainly through interactions.
- Metamodel-assisted SA. When each run is expensive, fit a surrogate (kriging, polynomial chaos, random forest) to a designed set of runs and compute the indices on the surrogate, reporting the surrogate's own validation error.
- Derivative-based measures and automatic differentiation for models where a gradient can be obtained cheaply.
| Method | Cost (model runs) | Use when |
|---|---|---|
| Local OAT | 2k for k factors | Quick check near a nominal design; nearly linear model |
| Regression / PRCC on LHS | Hundreds | Monotonic model; want a ranked list quickly |
| Morris screening | r(k+1), r ≈ 10–50 | Many factors, need to eliminate the unimportant ones cheaply |
| Sobol indices | N(k+2), N ≈ 103–104 | Non-linear model with interactions; a defensible, quantitative attribution is needed |
| Metamodel-based | 102 runs + surrogate | Each run costs hours; accept surrogate error in exchange |
For M/M/1, W = 1/(μ−λ), so ∂W/∂λ = 1/(μ−λ)2. At λ = 12, μ = 15 the elasticity is (λ/W)∂W/∂λ = λ/(μ−λ) = 4: a 1% rise in arrival rate lengthens the mean sojourn time by 4%. At λ = 14.5 the elasticity is 29. The sensitivity itself depends violently on the operating point — which is exactly why a single local sensitivity number, quoted without its nominal point, is meaningless.
5.4 Reporting sensitivity results
- State the nominal values and the assumed range or distribution of every input.
- State the output measure and the sampling design used.
- Give a tornado diagram (local) or a bar chart of Si and STi (global).
- Separate parametric from structural uncertainty — changing a distributional assumption or a model mechanism is not a parameter sweep, and should be reported as an alternative scenario.
- Conclude with the decision-relevant statement: which factors must be measured better, which can be ignored, and whether the recommendation is robust.
6. Unit Summary
- Distribution choice is a modelling decision, not a formality; never use an unbounded symmetric distribution for a duration.
- Kendall's notation names a queue; Little's law L = λW holds for essentially every stable system and is the first check on simulator output.
- M/M/1: ρ = λ/μ, L = ρ/(1−ρ), W = 1/(μ−λ). Queues explode non-linearly as ρ → 1.
- Pollaczek–Khinchine shows congestion depends on service variability (CV), not on load alone.
- Noise is characterised by its spectrum (white, pink, brown, shot); SDE noise scales with √Δt, and coloured noise comes from filtering white noise.
- PRNGs must have long period, good equidistribution, reproducibility and independent substreams; test with TestU01, and always record the seed.
- Variate generation: inverse transform (derive it!), acceptance–rejection with expected cost c, convolution, composition, Box–Muller, thinning for non-homogeneous Poisson, Cholesky for correlation.
- Output analysis: discard the warm-up, use independent replications or batch means, report confidence intervals, and compare systems using the interval on the difference with common random numbers.
- Sensitivity analysis: local elasticities are cheap but blind to interactions; Morris screens, Sobol indices attribute variance (Si vs STi ), metamodels make it affordable.
6.1 Formula sheet
| Quantity | Formula |
|---|---|
| Little's law | L = λW; Lq = λWq |
| M/M/1 utilisation, queue, time | ρ = λ/μ; L = ρ/(1−ρ); Lq = ρ2/(1−ρ); W = 1/(μ−λ); Wq = ρ/(μ−λ) |
| Pollaczek–Khinchine (M/G/1) | Lq = ρ2(1+CV2) / (2(1−ρ)) |
| Inverse transform | X = F−1(U); exponential: X = −(1/λ)ln U |
| Acceptance–rejection | accept if U ≤ f(Y)/(cg(Y )); expected trials = c |
| Box–Muller | Z = √(−2 ln U1) cos(2πU2) |
| Confidence interval | X̄ ± tn−1,1−α/2 s/√n |
| Replications for target half-width | n ≈ n0(h0/h*) 2 |
| Elasticity | (∂y/∂x)(x/y) |
| Sobol indices | Si = Var(E[Y|xi ])/Var(Y); STi includes interactions |
6.2 Practice questions
Short answer (2–3 marks each)
- Explain Kendall's notation with the example M/D/2/10.
- State Little's law and one condition under which it holds.
- Derive the inverse-transform generator for the exponential distribution.
- What is the expected number of trials per variate in acceptance–rejection, and why?
- Why must a Normal distribution not be used for service times?
- Distinguish local and global sensitivity analysis in one sentence each.
- What does a total Sobol index much larger than the first-order index indicate?
Medium answer (5 marks each)
- Jobs arrive at a server at 20/hour with exponential inter-arrival times; mean service time is 2.5 minutes, exponentially distributed. Find ρ, L, Lq, W and Wq, and verify Little's law.
- Explain the properties required of a pseudo-random number generator and describe one test used to validate one.
- Describe the acceptance–rejection method and apply it to f(x) = 2x on [0,1].
- Explain warm-up removal and the batch-means method for steady-state output analysis.
- Describe the Morris screening method and state what μ* and σ tell you.
Long answer (10 marks each)
- Discuss random variate generation in full: inverse transform (with proof), acceptance–rejection (with efficiency), convolution, composition and special methods, giving one derived example each.
- Explain queueing systems, Kendall's notation, Little's law and the M/M/1 results; then explain, using the Pollaczek–Khinchine formula, why variability rather than load alone drives congestion, and why simulation is needed beyond these models.
- Explain sensitivity analysis: purposes, local versus global methods, Sobol variance decomposition, metamodel-assisted analysis, and how results should be reported to a decision maker.
- Describe the complete statistical procedure for a stochastic simulation study, from input distribution fitting to reporting a defensible comparison of two system designs.
6.3 Further reading
- A. M. Law, Simulation Modeling and Analysis — chapters on input modelling, random number generation and output analysis.
- S. Asmussen and P. W. Glynn, Stochastic Simulation: Algorithms and Analysis, Springer.
- L. Devroye, Non-Uniform Random Variate Generation (freely available) — the reference for Section 4.
- A. Saltelli et al., Global Sensitivity Analysis: The Primer — for Section 5.
- D. Gross and C. M. Harris, Fundamentals of Queueing Theory — for Section 2.